\ 


* f 


X-590-73-190 

prsprint '' 


m <t:3 l 



A SEMI-ANALYTICAL METHOD OF 
COMPUTATION OF OCEANIC TIDAL 
PERTURBATIONS IN THE MOTION OF 

ARTIFICIAL SATELLITES 


(NASA-TM-X-66298) A SEMI- ANALYTICAL 
METHOD OF COMPUTATION OF OCEANIC TIDAL 
PERTURBATIONS IN THE MOTION OF ARTIFICIAL 
SATELLITES (NASA) -2ir p HC $^2v§ 

£ ( CSCL 22C G3/31 


N 73- 277 6 


5 


Unclas 
09912 J 


- / 




PETER MU SEN 


JULY 1973 






fvo 


x y 




3 


w 


;vN. 

,aV' 

V 



. .. 

GODDARD SPACE FLIGHT CENTER 

BREENBELT, MARYLAND 


N 


Reproduced by 

NATIONAL TECHNICAL 
INFORMATION SERVICE 

US Department of Commerce 
Springfield, VA. 22151 


U' 


/ 


A SEMI-ANALYTICAL METHOD OF COMPUTATION 
OF OCEANIC TIDAL PERTURBATIONS IN THE MOTION 
OF ARTIFICIAL SATELLITES 


Peter Musen 


July 1973 


Goddard Space Flight Center 
Greenbelt, Maryland 



BASIC NOTATIONS 


i, j, k - the basic unit vectors in the equatorial system, the x-axis is 
directed toward the vernal equinox and z-axis coincides with 
the axis of rotation of the Earth 

G - the gravitational constant 

M - the mass of the Earth 

6 - the sidereal time at Greenwich 

r' - the mean geocentric position vector of a water particle 
r - the geocentric position vector of the satellite 

A = |r - r'| - the distance between the satellite and the mean position of the 

particle of water 

r° - the unit vector in the direction of r 

v - the equatorial components of r° in the equatorial 
system 

v'. <*>' - the east longitude and the latitude of the water particle in the 
state of equilibrium 

a. - the right ascension of the satellite 

S - the declination of the satellite 

x - the density of oceanic water 

k 0 - the mean density of the Earth 

R - the radius of the mean level surface of the sea water 
R£(y', 0', t) - the height of water tide 

- Love numbers associated with the effects of loading of the 
mass of water 

// 



a - the seini -major axis of the orbit of the satellite 


„ - the mean motion of the satellite, n 2 a 3 = GM 
a, e, a. rr, i - the mean elliptic elements of the satellite 

P = R/a - the parallactic factor of the satellite 
1 - the mean anomaly of the satellite 

f - the true anomaly of the satellite. 
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A SEMI-ANALYTICAL METHOD OF COMPUTATION 
OF OCEANIC TIDAL PERTURBATIONS IN THE MOTION 
OF ARTIFICIAL SATELLITES 


INTRODUCTION 

The problem of tidal influence on the motion of artificial satellites has 
captured the imagination of theoreticians and observers. 

This interest was stimulated by the fact that from the observed magnitude 
of the tidal effects on the motion of satellites one can determine the elastic 
response of the earth as it is "seen" by the satellite from outer space. 

The pioneering theoretical and computational work in this domain was done 
by Kozai (1965), Newton (1968) and Kaula (1969). The problem is comparatively 
simple, theoretically and numerically, if only the influence of the "solid Earth" 
tides on the motion of the satellite is to be taken into consideration. A number 
of authors (Smith et al., 1971), (Anderle, 1971), (Douglas et al., 1972) have 
treated this problem with success. 

The electronic computer was applied recently to develop a semi-analytical 
theory of the Solid Earth tidal effects (Musen and Estes, 1972), (Musen and 
Felsentreger, 1973), in the form of Fourier series with purely numerical 
coefficients. However, recent observations indicate a disagreement with the 
theory and show an apparent dependence of elastic parameters of the Earth 
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(Love numbers) upon the orbital inclination of the satellite (Smith, 1973), 

(Douglas et al., 1972). The influence of the tides of the hydrosphere can be 
the primary cause of this disagreement and of the apparent dependence of Love 
numbers upon the inclination (Lambeck and Cazenave, 1973). It seems that the 
perturbations caused in the motion of satellites by the ocean tides can reach as 
much as 20% of the perturbations caused by the tides of the Solid Earth, and 
both kinds of perturbations have the same frequency spectrum. The most sig- 
nificant perturbations, which will affect the values of Love numbers, are caused 
by the ocean diurnal tides K 19 Oj , P a and by semi-diurnal tides M 2 , S 2 and P 2 . 
The Kj-tide contributes to the largest periodic term in the expansion of per- 
turbations in the orbital inclination and thus influences the value of Love number 
k 2 . In addition, the influence of normal modes of particular oceans also deserves 
consideration. The central computational problem consists of the determination 
of amplitudes and lags of the satellite's perturbations from the amplitudes and 
lags of the tidal constituents as exist over the surface of the world ocean. 

In the present article we discuss the method of expansion of the satellite's 
perturbations, as caused by the oceanic tides, into Fourier series with the 
arguments SI and rr of the satellite and£ ,V, F, D and r of the Moon. The coef- 
ficients in this expansion are purely numerical and peculiar to each particular 
satellite. Such a method is termed as semi- analytical in Celestial Mechanics. 
Gaussian form of the differential equations for variation of elements, with the 
right hand sides averaged over the orbit of the satellite, is convenient to use 
with the semi-analytical expansion. 
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EXPANSION OF THE DISTURBING FUNCTION 


The disturbing function associated with the direct attraction of the satellite 
by the tidal shell of water has the form: 

U = G/<R 3 JJ t) da< (i) 

where 

dcr' = cos 0'd0' dy', 

i = 0(1(T 7 ). 

The effects of loading on the bottom of the sea (Hendershott, 1972) are tempo- 
rarily omitted and will be introduced at the next stage. The double integration 
in (1) is performed over the mean surface of the global ocean (or over the mean 
surface of a particular sea whose tidal influence on the satellite we want to 
determine). 

Introducing the mean density of the Earth 


„ = 3M 
0 477 R 3 

we re-write (1) in the form: 


U = 


3 gm 

R 


K 


1_ 

477 


JJ ! 


t) do-' 


( 2 ) 


Substituting the expansion 


R 

A 


2 --err 

n = 0 m=0 


N nm P „m C s 1 n S) P nm ( s i n 0') cos m(a - 9 - y') 
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for m / 0 


N = + 1, 


N 


. (n - m) 
(n + m) 


into (2) and introducing Love numbers to account for the effects of loading 
(Lamb, 1945), we obtain for the oceanic tidal disturbing function: 


V 


GM 

IT 


+ 00 


c 1 + 


n = 0 



(3) 



N P (sin §) • 

nm nm v 7 


l 

4 TT 





t) P nm (sin 0') cos m(a - 6 - y') dcr' 


m = 0 S 

We assume that the tidal oscillations of the sea at a given point (y' , 0') of the 
mean sea surface can be represented as a sum of periodic constituents of the 
form: 


Z(y', 0', t) = Z(y\ 0') COS [(a)t + V) - S(y', 0')] 

i = 2z, (4) 

the phase angle s(y', 0') is the retardation of the maximum ("high water") of 
z(y', 0', t)at (y', 0') relative to its high water at Greenwich. The basic data 
plotted on the geographic maps of a constituent are the isochrones of its retarda- 
tion (the co-tidal lines) and the lines of equal amplitudes (the co-range lines). 


From these data one can compute the lag and the amplitude of the corre- 
sponding tidal term in the perturbations of the satellite. To simplify the 
formulas we represent a constituent not in the standard form (4) but rather in an 
"old fashioned" manner as: 

z (y'> 0'- t) = Z(y\ 0') cos [m (6 + y') + q(y', 0') + 0] • 
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where 

0 = crt + /? 

and o' is small. The local sidereal time <9 +y' disappears from the disturbing 
function in the process of elimination of the short period terms; the slowly 
changing part 0 + q remains. 

The most significant long period effects are produced in the motion of 
satellites by diurnal (m = 1) and semi-diurnal oceanic tides (m = 2). Among 
the diurnal constituents m Kj and S K 1 (combined together into ), O and are 
the most significant ones. Kj produces in the perturbations of the satellite 
terms with the argument <fl . Such a term has the largest amplitude in the tidal 
perturbations of i as caused by the solid Earth and it is proportional to Love 
number k 2 . The contribution by K x is of special significance in the case of 
satellites with very high inclination (Isis type). It is clear, therefore, that the 
oceanic tides with the period of one sidereal day affect considerably the deter- 
mination of k 2 . Two other constituents, Oj and P 1 , also have their share in 
modifying the value of k 2 . However, their contributions will have small ampli- 
tudes and shorter periods than the one previously mentioned. The arguments \p 
associated with O x and P : are 2V + 2D + 2T and 21' + 2r, respectively. The 
computation of the influence of the oceanic K ^constituent is facilitated by the 
presence of large amphidromies in the North Atlantic and the Pacific. The 
Ki-tides are highest in the Pacific and the Indian Ocean (Defant, 1961). 
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In order to fully understand the structure of the <0. -term in the tidal per- 
turbations of the satellite the numerical integration of the Laplace tidal differ- 
ential equations shall be repeated for Kj over the World Ocean to check or to 
improve Dietrich's results (1944). This integration shall be performed in a 
manner similar to one used recently for M 2 and S 2 , taking coastal boundary 
conditions and friction into consideration. Besides the real tides, the oceans 
are also subjected to free oscillations. For example, the North Atlantic has 
free oscillations with the periods of 21 h .2, 14 l ?0 and 11^5, with one, two and three 
amphidromic systems, respectively (Platzman, 1972). Oscillations of the sea 
such as these, with the periods of approximately a day or half-a-day, can also 
produce long period perturbations in the motion of the satellite. 


We substitute (5) into (3) and omit from the result all of the terms which 
contain local sidereal time in the argument. The resulting disturbing function 
is free from the short period perturbations caused by the Earth rotation 
(Lambeck and Cazenave, 1973). We have for a given constituent: 


where 


V 

m 


3. GM 

2 IT 


K 


+ oo 



n - m 


V - t~» n+ i 
nm 


P”*(l tk')AN _ 

v n ' n nm \ 


n+ 1 


P nm( sin S ) COS ( ma + '/ ;n )> 


(6) 


(7) 


- at + (/? + e n ), 
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and 


A n COS e n = ^ JJ z (y'. <£') p J sin <*>') COS q(v'. 0') dcr ' 

s 

A sine = — 

n n 477 

s 

n = m, m + 1 

As pointed out already we suggest only the computation of the influence of 
diurnal (m = 1) and of semi-diurnal tides (m = 2). The amplitudes A n and the 
lags e n are different for each constituent. In order to simplify the notations 
and to avoid the use of multiple indices we use in (8)— (9) the same notations for 
a given n, independently of m and of the constituent. The formulas (8)-(9) per- 
mit one to evaluate the amplitudes and lags of an oceanic tidal component in the 
perturbations of the satellite from the amplitudes and lags distributed over the 
surface of the ocean. The existence of amphidromic systems facilitates greatly 
this evaluation. The necessary basic information can be taken from the charts 
of the co-tidal and co-range lines. Such charts, obtained recently by integrating 

numerically Laplace tidal differential equations over the global ocean, are given 
for M 2 and S 2 components (without the effect of friction) by Bogdanov and Magarik 

(1967), for M 2 , including friction (mainly along the costs) by Pekeris and Accad 
(1969) and, most recently by Hendershott (1972). For the chart was given by 
Dietrich (1944). 


Z(y\ 4 >‘ ) P (sin $•) sin q(y\ <£') do-' 


(9) 
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We can express the spherical harmonics 


Y nra = P nm ( s 1 n &) COs(m, a + i/» n ) 


in (7) in terms of the components X. , /j. ,v of the geocentric unit vector of the 
satellite by using the general formula 


e ima P (sin S) 

nm x ' 


(n + m) ! 
2?m! (n - m) ! 


[( n-m)/2] 


(X + a. (X 2 + ^2)2j v n-m-2j 

j =0 


where 


( 10 ) 


i +i 


_ _ (n - m - 2j) (n - m - 2j - 1) 
4* (j + l)(m + j + 1) 


a . 
) 


(id 


= (- l) i+1 ( n ~ m K n - m - !)••• (n - m - 2j - 1) 

2 2j +2 . ( j + 1) ! (m + l)(m + 2 ) . . . (m + j + 1 ) 


j =0, 1,... 



or by means of the recursive relation 


(2n + 1) vY = (n - m + 


1 ) Y n + l,.+ 


(n + m)(X 2 + fj. 2 + v 2 ) Y 


( 12 ) 


n- 1 , m ' 


Such expansions in terms of X , fj. and v are convenient for the computation of the 
tidal perturbations in rectangular coordinates by means of numerical step by step 
integration. 
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However, the physical characteristics of the tidal perturbations, the relative 
importance of amplitudes and periods , as well as the presence of resonances , can 
be read and traced easier from a periodic expansion than from a numerical inte- 
gration. For this reason we develop in the present work the theory of the oceanic 
tidal perturbations based on the expansion into trigonometric series with the 
arguments tt and si of the satellite and -t, V , F, D and F of the Moon. The coef- 
ficients of the expansion are purely numerical. We make use of the Gaussian 
form of the differential equations for variation of elements so that the numerical 
values of the orbital inclination and eccentricity can be substituted from the outset. 
Such an approach to tidal perturbations, termed as "semi -analytical" in Celestial 
Mechanics, leads to a compact form of trigonometrical series with arguments 
linear in time. Operations with them can be handled efficiently on an electronic 
machine using the program developed at Goddard Space Flight Center by R. Estes. 

Substituting the value of i into 


\ = + cos 2 i c o s ( f + 77 ) + sin 2 - cos(f + tt - 2 Si), 
2 2 


M = + COS 2 i sin(f +7 7) - sin 2 - sin(f + 77 - 2 SI), 
^ 2 


(13) 


v - sin isin(f + 77-^1). 

We compute the semi -analytical expansions of 

Yjj = \ cos 0 - fi sin 0 


Y 91 = 3\v cos 0 - 3/U.v sin 0 
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and of 

Y 22 = 3(A 2 - /x 2 ) cos \p - 6 A/x sin 0 

Y 32 = 15v(A 2 - /x 2 ) cos i // - 30kfj.v sin 0 (15) 

into f ini te Fourier series with the arguments f,w , and 0. Making use of (12) 
we can then compute step by step the semi-analytical Fourier expansions of 

Y , y , . . . 

31’ Ml’ 

and 

y y 

M2’ x 52’ 

with the same arguments as before and numerical coefficients peculiar to each 
particular satellite. 


We can set 

A. 2 + M 2 + v 2 = 1 

in (12) when we compute Y nm , but we must keep this factor in the process of 
developing the recursive formula for the derivatives and for the force compo- 
nent normal to the orbital plane. To compute the general tidal perturbations 
normal to the orbital plane of the satellite we need the Fourier expansion of 


Z = R • V. Y 

nm u nm 


where V Q is the gradient operator with respect to r ° > 


V~J • B # a B 

v « = 1 a + J 5: + k s- 


and 


R • V. = + sin i sin ft sin i cosfl — + cos i — 

0 BA B/x 3v 


( 16 ) 
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By applying the V 0 -operator to (12) and setting X. 2 + fj? + v 2 =1 after the dif- 
ferentiation, we obtain: 


(2n + 1) zA7 0 Y nm + (2n + 1) kY nm 

= (n - m) V 0 Y n+l m + 2(n + m) r°Y n _ lm 

+ (n + m) V Y . 

' 0 n~ 1 , m 

Multiplying both sides of the last equation by R and taking 

k • R = cos i, r° • 'R = 0 
into account we deduce the recursive relation 

(2n + 1) vZ nm + (2n + 1) Y nm cos i 

= ( n - m > Z n+ l, m + ( n + m ) Z n-l, m (17) 

Making use of (14)-(16) we obtain 


Z n = sin i sin (fl +0) 

Z 2 1 = ^Zj! + 3Y U cos i ( 18 ) 

and 

Z 22 = 6\ sin i sin(f), +'/') + sin i cos(«n, + i//) 

Z 32 " 5vZ 22 + 5Y 22 COS 1 (19) 

Substituting the value of i for a given satellite and making use of (13) we expand 
(18)-(19) on the computer into Fourier series and making use of (17) obtain step 
by step the semi-analytical expansions of 


and 


Z 7 
31’ ^41’"' 


Z 7 

42’ ^5 2’ ' ‘ ' 
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After all necessary expansions are completed ip must be replacedby the proper value 
of <P . The formal accuracy of A. , p and v , with which we start the computations, 
is 10 7 - 10 8 . Because of the smallness of A n the final accuracy we need in the 
expansion of Y andZ nm , is only about 10 -3 - 10~ 4 . This indicates that the re- 
cursive relations for Y and Z serve our needs quite well, without a drastic 

nm nm 

loss of accuracy for a considerable number of steps. 

Unfortunately, we cannot say yet with complete security how many terms 
we require in the expansion of the disturbing function (7). The computations for 
particular satellites and their comparison with observations will show this more 
clearly. In connection with this point we wish to emphasize again the advantage 
of a semi-analytical way of treating the oceanic tidal perturbations of the satel- 
lite. The addition of any new term to the disturbing function can be done auto- 
matically and the decision about the necessary accuracy is left totally to the 
machine. 


DIFFERENTIAL EQUATIONS FOR VARIATION OF ELEMENTS 

At the present time we are able to observe only the long period oceanic 
tidal perturbations in the elements of the satellite. In the present exposition 
we make use of the Gaussian form of the differential equations for perturbations 
in elements and eliminate the short period perturbations by averaging the right 
hand sides of the equation over the orbit of the satellite: 


dSe 

UT 


« a 2 e 277 J 0 


(20) 
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d§77 /l - e 2 _ 1 

dt _,2 0 277 


na'e 


I 


27 r 


BH , /a 1 \3fl . J 
- a — — cos f + — + — sin f > dl 

. 3r r i _ e 2y Bf j 


n . 2 i dft 
+ 2 sin - , 

2 dt 


( 21 ) 


. . dSA 

s in l = + 

dt 


i i r 2n 

— • r(R • VQ) sin(f + 77 -XI) dl, (22) 

« aVl - e 2 277 J 0 


dS: 

"dt 


J C 2 77 

r (R* -VQ) cos ( f + 77 - A) dl 

o 


(23) 


dSL _ __2_ , J_ 
dt n a 2 277 


>277 


^0 


r — dl + (1 - /l - e 2 ) — + 2 /l - e 2 sin 2 - — , (24) 

9r 7 dt 2 dt 


As a result of eliminating the short period terms the semi-major axis a is not 
affected by tides. It is sufficient to consider the typical form of the disturbing 
function: 


Q = 3 GM _ K— y 

2 R nm 


or 


H = n 2 a 2 B _ Y 


n+ 1 


(25) 


where 


B = p n (1 + k') A N 

n m * x n ' r» 


n' n nm 
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From (25) we have 


rVfl = « 2 a 2 B nm ( — j [_ (n + 1) r° Y nm + rV r <> • V 0 Yj . 


Taking 


Vr° = A (I - r°r°) and r°-V.Y = nY 


0 nm nm 


into account we obtain: 


rVfi = n!a2B n.( j) ' - (2n + 1) r»Yj 


and, as a consequence, 


rR*-Vfi= »2 a 2 B (ff z . 


Substituting (25) and (26) into the differential equations (20)-(24) we reduce them 
to a form convenient for programming: 


d § e _ n /l - e 2 n 1 
, B • • - — 


ay* 1 

r) -A 


e / 1 - e 2 


ifC)' 


W dl 

nm 


0 . , i dSfl 

+ 2 sin^ - , 

2 dt 


W nm = (n + 1) MY + N — — 

nm 7 nm ^ £ 
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and 


M = i e + cos f + - e cos 2 f , 

2 2 


N = 2 sin f + - e sin 2f, 

2 


dSft 

~dt 


/ 1 - e 2 sin i ^ 7T J 0 


B ~ | (-] Z sin(f +Tr-fl)dl 


dSi n R 1 C /a \ n+1 7 m ,, 

dt nm 277 J Q [rj ™» 


(30) 


^=2nB (n + 1) • — f (-) 1 Y dl 
dt nmV ' 2tt J o \tJ nm 


(31) 


1 + /I - e 2 dt 


d §77 r 5 • 7 1 dS.fi> 

e + 2/1 - e sir - 

2 dt 


(32) 


The factors multiplying (a/r) n+1 in the integrands are finite Fourier series 
with the arguments f +tt, SI and y/ . Thus, the typical integrals to be evaluated 
have the form: 


^ ( 7 ) sin [k < f +7r > + + ^ dl 


X 0 n 1,k (e)- ^ (k77 + jft + \p). 
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where X~ n 1,k are Hansen coefficients (Tisserand, Vol. I). To compute these 
coefficients we can use the relation 

x r 1,k ( e ) = ( n k 1 )~(|) k (l - e2 ) _n+1/2 ( ss ) 


n +l k-n+2 
"2 2 


k + 1, e : 


For satellites presently used for the evaluation of tidal effects, we have e < 0.07 
and, consequently, we need only few terms in the hypergeometric series (33). 

If the eccentricity is very small then the presence of e as a "small divisor" in 
the right hand side of the differential equations might cause difficulty in expand- 
ing the semi-analytical theory of perturbations in elliptical elements. It is 
evident from (33) that such a difficulty can arise only if k = 0, i.e., if there are 
terms of the form 



cos 

sin 


(ja + </o 


under the integral sign. The terms of this type disappear in (27) because of 
differentiation with respect to f. As a result a semi-analytical expansion of the 
tidal dSe/dt can always be arranged and in this case e is only a spurious 
"small divisor." It does not cause any real difficulty. The real difficulty 
associated with small e can arise only on developing the tidal dS^/dt. We can 
see from the form of M, N and W nm that e will appear as a real "small divisor" 
only if there is a term of the form 
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cos . 

K sin ( f + ir + ]il + <A) 

in the expansion of Y nm . This will happen only if n is an odd integer. For n 
even the semi-analytical expansion of the tidal dSw/dt always can be arranged. 

The arguments in the resulting series for derivatives of the elements are all 
linear in time and, consequently, the integration is straightforward (except in a 
case of a sharp resonance). 

CONCLUSION 

In the present work we give a set of formulas necessary to expand the satel- 
lite's oceanic tidal perturbations in a semi-analytic form, as Fourier series 
with purely numerical coefficients, using a computer. All arguments are linear 
with respect to time (providing we do not have sharp resonances). This form 
permits one to compute the tidal perturbations and at the same time to make a 
judgement about the relative importance of amplitudes and periods of different 
terms. We select a semi-analytical, and not an analytical method, because of 
its compactness and because the basic information about the oceanic tides is 
given either in the semi-analytical form or in the form of charts of co-tidal 
and co-range lines. 

On the basis of already existing information we can obtain a general estimate 
of the oceanic effects on the motion of an artificial satellite. 

However, in order to fully understand these effects and to have a more accurate 
information about them we have to obtain more complete information about the 


global oceanic tides themselves. 
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We suggest a new numerical integration of Laplace tidal differential equa- 
tions for those tidal constituents for which the work was done either on a limited 
basis or long ago. 

In particular, it will be of great interest to obtain new information for K,. 
This constituent contributes to the biggest amplitude in the expansion of the 
tidal perturbations of the orbital inclination. The corresponding term has the 
right ascension of the node as the argument. 

The knowledge of oceanic perturbations in the motion of satellites will per- 
mit us to improve the values of the elastic parameters of the earth. 
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